Cumulant expansion for fast estimate of non-Condon effects in vibronic transition 

profiles 
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When existing, cumulants can provide valuable information about a given distribution and can 
in principle be used to either fully reconstruct or approximate the parent distribution function. 
A previously reported cumulant expansion approach for Franck-Condon profiles [Faraday Discuss., 
150, 363 (2011)] is extended to describe also the profiles of vibronic transitions that are weakly 
allowed or forbidden in the Franck-Condon approximation (non-Condon profiles). In the harmonic 
approximation the cumulants of the vibronic spectral profile can be evaluated analytically and 
numerically with a coherent state-based generating function that accounts for the Duschinsky effect. 
As illustration, the one-photon 1 ^Ag — >■ 1 ^B2u UV absorption spectrum of benzene in the electric 
dipole and (linear) Herzberg- Teller approximation is presented herein for zero Kelvin and finite 
temperatures. 
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I. INTRODUCTION 

Vibrationally resolved electronic spectra {e.g. one- 
photon absorption and emission spectra) are within 
the Born-Oppenheimer framework usually interpreted in 
terms of Franck-Condon (FC) factors (FCFs) Ac- 
cordingly, one can try to obtain the shape of the spectral 
profile for a FC-allowcd transition from computed FCFs 
in frequency domain. However, the evaluation of FCFs 
for large molecular systems is challenging even within the 
harmonic approximation if one has to take Duschinsky 
mode mixing (rotation) into account. This is because 
multi-variate Hermite polynomials have then to be eval- 
uated for each FC integral, rather than only uni-variate 
Hermite polynomials as is the case for the comparatively 
simple parallel harmonic oscillator model. The computa- 
tional task becomes more difficult as the molecular size 
and temperature increases because the number of FC in- 
tegrals grows vastly. To describe FC-forbidden or weakly 
allowed transitions, one has to go beyond the Condon ap- 
proximation and employ for instance a Herzberg- Teller 
(HT) expansion Q of the electronic transition moment 
with respect to the normal coordinates. As a result, the 
calculation of the vibronic spectrum for a non- Condon 
process is even more difficult than for a FC-allowed tran- 
sition because one has to evaluate many matrix elements 
of the non-Condon operators which require for each HT 
integral in general the calculation of combinations of sev- 
eral FC integrals. 

The number of FC integrals and matrix elements of 
non-Condon operators to be evaluated in a sum-over- 
states approach can be significantly reduced with the help 
of rigorous integral prescreening strategies 0-0] ■ How- 
ever, this time-independent (TI) calculation of the spec- 
tral profile in frequency domain is still considerably more 



expensive than an alternative time-dependent (TD) ap- 
proach that exploits time-correlation functions (TCFs) 
(see e.g. Ref. Q), but offers the ability to directly assign 
individual peaks in the spectrum. As we have outlined 
earlier [1, 0, , a unified coherent state-based generat- 
ing function (CSGF) approach [l3| can be used both for 
rigorous integral prescreening strategies and TCF calcu- 
lations which combine the strengths of both approaches 
and complement each other favourably. Even in the less 
demanding TD approach, however, one usually has to in- 
vest significant computational time for often unnecessary 
spectral details. 

Cumulants (or moments) of a distribution (see e.g. 

Rcfs. [iM3) 

can deliver highly useful information. 
From this one can either attempt to reconstruct the spec- 
tral shape or try to estimate the relevant spectral profile, 
which can be exploited in subsequent TI and TD ap- 
proaches [l^. Cumulants of the vibronic spectrum can 
be obtained from the CSGF directly without computing 
the total spectrum in frequency domain. This method 
was exploited already in Ref. [13 for FC-allowcd transi- 
tions, and we report herein an extension of this method 
to incorporate non-Condon transitions. As application, 
we present the profile of the ^Aig — !■ ^B2u transition 
of benzene, which is in the electric dipole approxima- 
tion Franck-Condon forbidden, at various temperatures 
within the linear HT approximation and compare the cu- 
mulant expansion result with the TCF approach. 



II. METHOD 

The spectral profile {g{huj]T)) can be expressed via 
the Fourier transform (FT) of the TCF (x(t; T)) that 
depends on the time t and temperature T, namely 
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where hu is the transition energy and hujQ the 0' — tran- 
sition energy. The corresponding occupancy representa- 
tion for the TCF can be obtained from Fermi's Golden 
Rule, i.e. 
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where we have assumed an electric dipole transition with 
the electronic transition dipole moment {jj,{Q)), which is 
a function of normal coordinates of the initial electronic 
state. The TV-dimensional harmonic oscillator eigenstates 
of the initial and final electronic state are denoted by \v) 
and \v') with the corresponding harmonic energy vectors 
e and e^, respectively. The vibrational Hamiltonians of 
the the initial and final electronic states are H and H' , 
respectively, is the Boltzmann constant. -Ee'.e is the 
vibronic transition energy with respect to the 0' — tran- 
sition energy. The time-dependent part with the vibronic 
transition energy corresponds to a Dirac delta distribu- 
tion in frequency domain and accordingly this expres- 
sion is usually evaluated in frequency domain, which is 
straightforward. The spatial representation of the TCF 
in closed form can be found by evaluating the following 
quantum mechanical traces 



Tr^^ACQ)* oxp(-iH't/?t)A(Q) cxp(i//t/'»)exp(-/f/(fcBT))j 



(3) 



The traces can be evaluated with any complete basis. In 
our work coherent states were used (see e.g. Refs. [H, 
0, Q ) with the Duschinsky relation between initial and 
final state normal coordinates {Q' = SQ + d where S and 
d are the Duschinsky rotation matrix and displacement 
vector, respectively and Q' are the normal coordinates of 
the final state) and the HT expansion of the electronic 
transition dipole moment (fi ~ fi{0) + iJ-'^Qk where 

/i^ is the first derivative of jl with respect to Qk-)- 

The vibronic spectral density function can be related 
to a probability density function (PDF) because the vi- 
bronic transition process follows a certain PDF. If all cu- 
mulants or moments of a PDF are defined and available, 
the PDF can be reconstructed as follows 



x(t,T) = £itotexp ^ 



Vfe=l 



kl 



(4) 



where {E^, ^)'^{T) is the fc-th order cumulant at temper 
ature T. The cumulants of the spectral density func 
tion are normalised to give the total intensity gtot = 

|2 t n'^ I ,,' |2 efc 



coth( 2fc^j ). Moments (cumulants 
and moments are inter-convertible) can be obtained by 
partial derivatives of x with respect to time. 



t=o 



(5) 



Thus, cumulants can be evaluated analytically or numeri- 
cally by evaluating partial derivatives of x in Eq. (U) with 
respect to the time variable at i = 0. Analytic evaluation 
of the cumulants to arbitrary order within the linear HT 
approximation can be performed along the lines of the de- 
velopment in Refs. for the cumulants of FC profiles 
to arbitrary order. For numerical evaluation of low-order 
cumulants one needs to compute x only at the first few 
time steps. To obtain the corresponding moments numer- 
ically, Re{x{t,T)) and Im(x(i,r)) are used for the cal- 
culation of even and odd moments, respectively, because 
Q-iE,,,j/h ^ cos{E,,,,t/h) -isui{E,,,,t/h) in Eq. ^ (see 
e.g. Ref. 0). 

The closed form of x(t, T) within the linear HT approx- 
imation can be found in Ref. [3, . Therein, additional 
flexibility was introduced to the GF by distinguishing 
individual vibrational modes to allow for rigorous pre- 
screening strategies and detailed control of the dynamics. 
This can be done by assigning different time and temper- 
ature variables to each vibrational degrees of freedom. 
The corresponding GF in an occupancy representation 
reads as follows 



G^(z,z';B,B')(/~'^)-AA|(0'|0)|-2 



9W 
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(6) 
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where the general operators / and g, which can be prod- 
ucts of momentum and position operators, are given in- 
stead of fi{Q). The time variables are assigned to the 
matrices z = diag(e''i*i/(2'^), . . . , e''"*"/^^'*') and z' = 
diag(e~"i*i/'^''\ . . . , e~'^«*"/(^''^) for initial and final vi- 
brational modes, respectively. Different temperatures 
can be given to the initial and final vibrational degrees of 
freedom via B = diag(l/(fcBTi), . . . ,1/(^3?^)) and B' 
= diag{l/{kBT{), . . . ,l/{kBT'^)), respectively. TV is the 
corresponding normalising factor. 

The electronic 1 ^Ag — I ^B2u transition of ben- 
zene is FC-forbidden in the electric dipole approxima- 
tion (/i(0) = 0) such that only the HT terms contribute 
to the spectral density function. The corresponding TCF 
is given as follows, here with same time and temperature 
(t,T) for all vibrational degrees of freedom. 



XFCHT(i; T) = |(0'|0) r E ■ h' G(A; T)'-^-^^ 



(7) 



III. RESULTS AND DISCUSSION 

The vibronic profiles for benzene's 1 ^Ag — > 1 ^B2u 
transition at zero Kelvin and finite temperatures are 
calculated with two methods, namely TCF and time- 
independent cumulant expansion (CE). We use herein 
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FIG. 1. Left part of the figure: The dashed Unes are drawn for 
the TCF-FFT approach with a Lorentzian fine shape function 
with FWHM of 50 cm^^. A time increment of 0.51 fs and 
a grid with 2^^ grid points are used for the corresponding 
FFT calculations. The experimental UV absorption spectrum 
of Fischer [2^ is additionally shown in read. Right part of 
the figure: The dashed lines are drawn for the TCF-FFT 
approach with a Gaussian line shape function of with FWHM 
of 500 cm~^. A time increment of 0.10 fs and a grid with 2^^ 
grid points are used for the corresponding FFT calculations. 
Solid lines are drawn for the curve obtained by Edgeworth 
expansion using up to 4-th order cumulants and, for 1000 K 
by Edgeworth expansion using up to 8-th order cumulants. 



the term time-independent CE, which was employed in 
R-cf. to distinguish this CE from the conventional 
(time-dependent) CE (see e.g. Refs. [U-[I3) which in- 
volves time integration for the cumulant calculation. The 
spectra obtained are displayed in Fig. [TJ To compute 
the vibronic spectra via the TCF method, the FFTW li- 
brary [U is used for fast Fourier Transform (FFT) . The 
approximate curves a re g enerated for the CE with Edge- 
worth expansion [l9l . [22| using the computed low-order 
cumulants. Some of the problems related to this type of 
expansion for the description of EC profiles are discussed 
in Ref. [l^. The moments (Eq. are calculated both 
analytically and numerically, the latter by taking partial 
derivatives of x in Eq. [3]with respect to time. Analytical 
and numerical results are compared in Table ID Required 
input data from electronic structure calculations for ben- 
zene, i.e. molecular equilibrium structures and corre- 
sponding harmonic force fields for each electronic state 
(^Ag and ^B2u) as well as first derivatives of the elec- 
tronic transition dipole moments arc taken from Ref. (23j 



(CASSCF/DZV). These data have recently been com- 
pared to results obtained via analytical derivatives tech- 
niques for electronic transition dipole moments within a 
time-dependent density functional theory framework [23l ■ 
The vibronic structure methods employed in the present 
work are implemented in a development version of our vi- 
bronic structure program package hotFCHT [sI-ItI. [isl [23| . 

The left hand side in Fig. [1] shows vibronic profiles 
from TCF-FFT which are convoluted by a Lorentzian 
line shape function with full width at half maximum 
(FWHM) of 50 cm~^ at temperatures elevating from K 
to 1000 K. This EC-forbidden vibronic transition is medi- 
ated by the non-totally symmetric vibrational modes in 
the irreducible representation e2g. The main feature of 
the vibronic spectrum is from progressions in the totally 
symmetric C-C stretching mode (963 cm~^) building on 
the so-called false origin from a single excitation of a 
non-totally symmetric (e2g) in-plane bending mode (575 
cm~^) as indicated in the spectrum at zero Kelvin. The 
calculated spectrum at 300 K is compared with the exper- 
imental data of Fischer [2^ . The two spectra agree fairly 
well in the low energy region but the computed peaks at 
higher energies are slightly shifted to larger wavenumbers 
due to the harmonic approximation. As temperature in- 
creases the vibrational structure becomes very congested 
and washed out. At 1000 K (only employed for testing 
the method) one can not see a resolved vibrational struc- 
ture any longer, only the corresponding envelope. 

On the right hand side of Fig. [T] the two methods 
(TCF-FFT and CE-Edgeworth) are compared for in- 
creasing temperatures. The spectra are convoluted in 
the TCF-FFT curves (dashed lines) with a Gaussian 
line shape function of 500 cm^^ for FWHM. The sec- 
ond moments (4.51x10^ cm^^ (/icq)^) of the Gaussian 
line shape function is added to the second moments of 
the vibronic spectrum to take the line shape function 
into account (see Ref. [l^ for the rationalisation and de- 
tails). The relatively broad line shape function is used 
for the TCF-FFT curves to have vibrationally relatively 
structureless spectra for comparison. At 0, 300 and 500 
K, the TCF-FFT curves still show vibrational structure 
and the CE-Edgeworth curves (solid lines) look like non- 
linear regression curves of the corresponding TCF-FFT 
spectra. When the vibrational structures are also essen- 
tially smoothed out in the TCF-FFT curves at 1000 K, 
the two approaches agree with each other extremely well. 
Up to the 4-th order cumulants arc used for 0, 300, 500 
K and up to 8-th order cumulants are computed for 1000 
K. 

In table U the moments computed numerically (via 
derivatives) and analytical are compared. At low orders 
and all temperatures the two methods agree well and 
for higher orders still the agreement is satisfactory. One 
of the advantages of the numerical method is that one 
only needs to compute the TCF for the first few time 
steps. The analytical method usually meets a combina- 
torial problem in high order cumulant calculations due to 
the analytic derivatives of the inverse matrix. The second 



4 



advantage is that it is easy to include linear and nonlinear 
non-Condon effects. The third advantage is that one can 
incorporate general line shape functions which would not 
have well defined cumulants (see the discussion on page 
415 of Ref. HI). 

TABLE I. Analytically and numerically computed cumulants. 
4.51 X 10* cm~^ (hco)^ is added to the second moments to take 
the Gaussian line shape function (FWHM = 500 cm^^) into 
account; see Ref. [l^ for details. A time increment of 0.10 fs 
is used for computing the numerical derivatives. 

(g.",.)/(<:m-' har 

T = OK r = 300 K T = 500 K T = 1000 K 

Analytical Numerical Analytical Numerical Analytical Numerical Analytical Numerical 

1 2.61x10^ 2.61x10' 2.47x10' 2.47x10' 2.12x10' 2.12x10' 1.12x10' 1.12x10' 

2 9.24x10" 9.21x10° 8.64x10° 8.62x10° 7.38x10° 7.36x10° 5.64x10° 5.63x10° 

3 4.07x10" 4.05x10'° 3.77x10'° 3.75x10'° 3.17x10'° 3.15x10'° 2.03x10'° 2.01x10'° 

4 2.14x10" 2.12x10" 1.97x10" 1.95x10" 1.66x10" 1.64x10" 1.29x10'' 1.27x10" 



5 - - - - - - 8.00x10" 7.82x10" 

6 - - - - - - 6.55x10^' 6.33x10^' 

7 - - - - - - 5.84x10^° 5.54x10^° 

8 - - - - - - 6.16x10^° 5.71x10^° 



TABLE IL Mean excitation wavenumbers of individual vibra- 
tional modes in e2g 



Mode (cm 


') K (cm" 


') 300 K (cm" 


') 500 K (cm" 


') 1000 K (cm"') 


57.5 


2.62x10^ 


2.90x10^ 


3.79x10^ 


7.25x10^ 


575 


2.62x10^ 


2.90x10^ 


3.79x10^ 


7.25x10^ 


1237 


1.60x10' 


1.78x10' 


4.64x10' 


2.51x10^ 


1237 


1.60x10' 


1.78x10' 


4.64x10' 


2.51x10^ 


1665 


2.54x10' 


2.38x10' 


3.09x10' 


1.65x10^ 


1665 


2.64x10' 


2.38x10' 


3.09x10' 


1.65x10^ 


3389 


8.14x10' 


7.49x10' 


6.08x10' 


6.37x10' 


3389 


8.14x10' 


7.49x10' 


6.08x10' 


6.37x10' 



Mean excitation wavenumbers are computed for the 
HT active e2g symmetric vibrational modes of the final 
state. The mean excitation energy can serve as a parame- 
ter for the individual vibrational degrees of freedom as an 
effective reorganisation energy or a Huang-Rhys factor 
(when normalised by its harmonic energy), which can be 
characterised as a function of structural deformation, fre- 
quency change, Duschinsky mode coupling and tempera- 
ture both in the Condon and non-Condon approximation. 
One might naively expect a larger mean excitation energy 
as the temperature increases but the mean values of high 
frequency modes (1665 and 3389 cm~^) in some interme- 
diate temperature rages arc smaller than at zero Kelvin. 
This can be rationalised as follows: Because the Duschin- 
sky mode mixing between low and high frequency modes 
is small in the present case, the high frequency modes can 
not obtain thermal energy from the low frequency modes 
efficiently. Thus the high frequency modes are almost 
thermally inactive, whereas the total intensity ((?tot) in- 



creases as temperature increases. In the mean energy 
calculation of the high frequency modes at finite tem- 
peratures the denominator (total intensity) increases be- 
cause low frequency modes accept thermal energy while 
the numerator (excitation of high frequency modes) stays 
constant. Therefore the mean excitation energies of high 
frequency modes are reduced at finite temperatures. If 
Duschinsky rotation couples the low and high frequency 
modes significantly (see e.g. Ref. (H), however, thermal 
energy can be transfered to the high frequency modes via 
the low frequency modes in the initial state, accordingly 
the mean excitation energies of high frequency modes can 
increase as temperature increases. 



IV. CONCLUSION AND OUTLOOK 

We have discussed a cumulant expansion method for 
describing non- Condon transitions and applied it to the 
prototypical one-photon electric dipolc 1 ^Ag — 5- 1 ^B2u 
transition of benzene, which is FC forbidden but HT 
allowed in the linear HT approximation. The method 
is particularly powerful when one does not require all 
the details of the vibronic structures, but rather only 
quantities such as peak maximum, mean and variance 
of the spectral shape. This method is computation- 
ally much cheaper than the sum-ovcr-states and time- 
correlation function approach. Moreover, the informa- 
tion (e.g. relevant energy window) from the cumulant 
expansion method can be used in the calculation within 
the other two methods. Herein we compared a numerical 
approach for the calculation of cumulants with the results 
from an analytical scheme. The results obtained numer- 
ically are still fairly good. With this method, one can 
incorporate easily nonlinear non-Condon terms and vari- 
ous line shape functions. In the time-correlation function 
calculation the real part and imaginary part at each time 
step provide automatically the even and odd moments, 
respectively. In the first few time steps we already have 
the first few moments available and the probability distri- 
bution function (information) becomes complete as time 
progresses. 
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